#!/usr/bin/python
# This program calculates average quality values of Qual files
# Can be single or multiple quality fasta files
# Usage: qualityStats.py filename.qual

import sys
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

qualfile = open(sys.argv[1], "rU")
MINQUAL = int(sys.argv[2])

avgqual = 0
qualrecs = 0

qlist = []

for q in SeqIO.parse(qualfile, "qual"):
    rec = q.letter_annotations["phred_quality"]
    total = 0
    for i in rec:
        total += i
        qlist.append(i)
    avgqual = float(total)/len(q)
    qualrecs += 1

print "There are " + str(qualrecs) + " quality files in " + sys.argv[1]
print "Average quality = ", avgqual

belowmins = 0

for i, j in enumerate(qlist):
    if j < MINQUAL:
        belowmins += 1

print "There are " + str(belowmins) + " bases below " + str(MINQUAL)

qualfile.close()

